clc; clear; close all;

% 参数设置
g11 = 1.4; g12 = -2; g21 = 1.4; g22 = 1.6;

% 参数网格
N = 1000;
eta_list = linspace(-3, 3, N);
k_list = linspace(-8, 8, N);

% 颜色映射，顺序和分类函数对应
abbrs = ["SNP", "UNP", "USP", "SNF", "UNF", "HBP", "FBP"];
colors = [
    0 0 0;       % SNP 黑
    0 1 0;       % UNP 绿
    0 0 1;       % USP 蓝
    0.5 0.25 0;  % SNF 棕
    1 0 0;       % UNF 红
    0 1 1;       % HBP 青
    1 1 0;       % FBP 黄
];

% 结果矩阵（字符串类型）
TYPE = strings(N, N);

% 并行计算稳定性类别
parfor i = 1:N
    eta = eta_list(i);
    local_type = strings(1, N);
    cos_pi_eta = cos(pi * eta);
    for j = 1:N
        k = k_list(j);
        % 计算 c1, c2
        c1 = 2 - g11 - g22 - k * cos_pi_eta;
        c2 = (1 - g22) * (1 - g11 - k * cos_pi_eta) - g12 * g21;
        
        % 分类
        local_type(j) = Classify(c1, c2);
    end
    TYPE(i, :) = local_type;
end

% 统计计数
counts = zeros(1, length(abbrs));
total_points = N*N;
for c = 1:length(abbrs)
    counts(c) = sum(TYPE(:) == abbrs(c));
end

% 输出占比
fprintf('稳定性类别占比统计（总点数 %d）：\n', total_points);
for c = 1:length(abbrs)
    percentage = counts(c) / total_points * 100;
    fprintf('%s: %d 个点，占比 %.2f%%\n', abbrs(c), counts(c), percentage);
end

% 绘图
figure;
hold on;
[X, Y] = meshgrid(eta_list, k_list); % 注意顺序，X对应eta，Y对应k

for c = 1:length(abbrs)
    mask = (TYPE' == abbrs(c));  % 转置对应meshgrid索引
    scatter(X(mask), Y(mask), 8, repmat(colors(c,:), sum(mask(:)), 1), 'filled');
end
xlabel('\eta', 'FontSize', 14);
ylabel('k', 'FontSize', 14);
grid on; box on;
legend(abbrs, 'Location', 'best', 'Orientation', 'horizontal', 'FontSize', 12);
hold off;

% 稳定性类别占比统计（总点数 1000000）：
% SNP: 91802 个点，占比 9.18%
% UNP: 149064 个点，占比 14.91%
% USP: 67676 个点，占比 6.77%
% SNF: 166126 个点，占比 16.61%
% UNF: 489562 个点，占比 48.96%
% HBP: 22072 个点，占比 2.21%
% FBP: 13698 个点，占比 1.37%